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Abstract: This paper investigates experimental means of measuring the 

transmission matrix (TM) of a highly scattering medium, with the simplest 
optical setup. Spatial light modulation is performed by a digital micromirror 
device (DMD), allowing high rates and high pixel counts but only binary 
amplitude modulation. We used intensity measurement only, thus avoiding 
the need for a reference beam. Therefore, the phase of the TM has to be 
estimated through signal processing techniques of phase retrieval. Here, 
we compare four different phase retrieval principles on noisy experimental 
data. We validate our estimations of the TM on three criteria : quality 
of prediction, distribution of singular values, and quality of focusing. 
Results indicate that Bayesian phase retrieval algorithms with variational 
approaches provide a good tradeoff between the computational complexity 
and the precision of the estimates. 
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1. Introduction 


Wave propagation in complex media is a fundamental problem in physics, be it in acoustics, 
optics, or electromagnetism HI. In optics, it is particularly relevant for imaging applications. 
Indeed, when light passes through a multiply scattering medium, such as a biological tissue 
or a layer of paint, ballistic light is rapidly attenuated, preventing conventional imaging tech¬ 
niques, and random scattering events generate a so-called speckle pattern that is usually con¬ 
sidered useless for imaging. Recently, wavefront shaping using spatial light modulators (SLM) 
has emerged as a unique tool to manipulate multiply scattered coherent light, for focusing or 
imaging in scattering media O. In essence, these methods use the linearity and time-reversal 
symmetry of the wave propagation, whatever the complexity of the medium, to control the out¬ 
put speckle field, by manipulating the light beam impinging on the scattering sample. Different 
wavefront shaping approaches rely on digital phase-conjugation nil or iterative algorithms 
(H, but it is also possible to measure the so-called transmission matrix (TM) of the medium ||6l, 
which fully describes light propagation through the linear medium, from the modulator device 
to the detector. This approach has been particularly efficient for focusing, imaging cmi and 
for studying the transmission modes of the medium O- These methods are not only valid for 
scattering material but can also be applied to other complex transmission system, most notably 
multimode fibers, turning them into minimal footprint endoscopes ifTOlfTTlfT^fTSl . 

A major limitation of most of these techniques for imaging is their speed. Indeed, the wave- 
front shaping process must be faster than the stability time of the medium, which can be of 
only a few milliseconds in biological tissues. Yet, most of the works reported so far have relied 
on phase modulators which are usually slow (few tens of Hertz for liquid crystal modulators). 
Micro Electro-Mechanical Systems (MEMS) modulators are much faster, but are usually not 
phase-only. As a promising alternative for wave shaping in complex media. Digital Micromir¬ 
ror Device (DMD) technology ina offers binary amplitude modulators {i.e., ON or OEE) op¬ 
erating at > 20kHz, with high pixel counts (10^) and low pitch (around 10 microns), all this at 
low cost. These binary amplitude modulators have been used as phase modulators, using ap¬ 
propriate diffraction and filtering, e.g. by Lee-type amplitude holography 03 na, as shown on 
Eig. [^). While phase control is more effective for wavefront shaping than amplitude control, 
some works reported on using DMD as genuine binary amplitude modulators for wavefront 
shaping through opaque scattering media, albeit usually yielding lower overall efficiency than 
phase modulators for focusing or mode matching EllIlllISl . The DMD configuration can also 
be optimized using genetic algorithms ll2Qll to maximize the intensity enhancement. 

Eor the measurement of a TM, an additional issue lies in accessing the amplitude and phase 
of the output field, that in optics usually requires a holographic measurement, i.e. a reference 
beam, as shown on Eig.[^). This reference beam can either be co-propagating in the medium 
EEIl, or use an external reference arm Eaiii. The phase and amplitude of the measured 
field can then be extracted by simple linear combinations of interference patterns with a phase- 
shifted or off-axis reference. This however poses the unavoidable experimental problem of the 
interferometric stability of the reference arm. 

In this work we report on the full measurement of the complex TM of a multiply scatte¬ 
ring medium, using a DMD binary amplitude modulator as an SLM, with no reference on the 
detection side, as shown on Eig. [T]:). This approach combines the high-speed and high pixel 
counts allowed by DMD devices, with the simplicity and robustness of a reference-less optical 
setup. However, it involves advanced signal processing algorithms for phase retrieval, run on a 
sufficiently large number of input-output calibration measurements. In this study, we compare 
the performance of four phase-retrieval algorithms |[23l|24l[25l[26l, for the estimation of a TM 
based on actual noisy experimental measurements. We assess their performance as a function 
of the number of measurements, and compare their relative computational cost. We then show 
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Fig. 1: Different experimental approaches for measuring the complex-valued transmission ma¬ 
trix of a scattering medium with a binary DMD amplitude modulator, (a) use of a reference 
arm for retrieving the phase of the output field by off-axis or phase-shifting holography; (b) 
Using the DMD as a spatial phase modulator by displaying amplitude holograms, and using the 
unmodulated parts of the field as a phase-stable reference; (c) the presented approach where 
only intensity values are measured. 


that the distribution of the singular values of the measured TM varies according to random ma¬ 
trix theory. Finally, we demonstrate that single- or multi-point light focusing can be achieved, 
using an -regularization algorithm (271 for determining the optimal DMD binary input pat¬ 
tern. In addition to being an interesting signal processing problem, our approach is particularly 
relevant for real-life applications of the TM approach, since it allows a simple, fast and robust 
implementation. 

2. Experimental setup 

Our experimental setup, described in Fig.[^ uses a DMD-array from Texas Instrument (1920 x 
1080 tilting micromirrors), driven by the DTP V-9500 VIS module (Vialux). The DMD is made 
of mirrors that can switch between two angular positions separated by 24°, thus reflecting each 
pixel either toward a beam dump (pixel OFF) or towards the focusing system (pixel ON). 

Under Matlab, an amplitude mask is computed and loaded on the DMD. The pattern corre¬ 
sponding to the ON pixels is focused on the surface of a thick scattering medium by means of 
a / = 100 mm lens LI (thus the DMD pixels correspond roughly to incidence angles on the 
sample). The sample is a ^ 100 microns thick layer of white paint, which is thick enough in 
order to considerably mix the light on the other side, producing a complex speckle interfering 
pattern. This speckle pattern is collected through a microscope objective (L2) and detected on 
a camera (AVT Pike F-IOOB). In order to measure the TM, we need to send a large series of 
input patterns (typically a few times the number of input pixels we wish to control), in a time 
over which the medium can be considered stationary. For this purpose, we use the “high speed” 

































































driver provided with the DMD in order to load all the to-be-projected random amplitude masks 
to the memory of the DMD driver module, and we trigger the display of each mask via a DAQ 
card (National Instruments, PCI-6221) and a waveform generator. In the same way, in order to 
be as fast as possible, we also only consider a subregion on the camera of size of 400 x 400 
pixels. The overall acquisition rate is 31 images/second. To monitor the stability of the medium, 
we periodically measure the correlation of the speckle image corresponding to the same input 
mask. We therefore quantify the stability of the medium, which is better than 98% over the total 
measurement time (typically around 5 minutes). 



Fig. 2: Experimental scheme: A 532 nm CW laser is expanded through a telescope in order 
to obtain an homogeneous beam. Through a rectangular mask, it illuminates the DMD which 
acts as binary amplitude spatial light modulator. The DMD reflects the light in two different 
directions corresponding to either ON (unit transmission) or OFF (the light is deviated towards 
a beam dump). The transmitted pattern is focused by a first lens LI on the scattering medium - 
here a white paint layer -, acting as a thick multiply scattering medium. The transmitted speckle 
pattern is collected by a microscope objective and is observed through a polarizer P on a CCD 
camera. 


3. Estimating the TM with intensity-only measurements and binary inputs 

The measurement of the TM can be formalized as a calibration problem: given P incoming 
waves, assumed perfectly known, which model explains at best the observed outputs? In our 
case, this inverse problem reduces to the well-known problem of phase retrieval. 

Let G {0,1}^ stand for the binary DMD inputs related to the /i-th acquisition, where N 
is the number of pixels (mirrors) used on the DMD. We assume that the partial observations of 
the sole moduli of the transmitted waves (the square root of the camera measured intensities), 
denoted by y^ G , obey 

yf, = |Dx^|, (1) 

where D is the TM complex-valued transmission matrix characterizing the scattering material, 
and M is the number of observed pixels on the camera. 






Then, adopting a matrix formulation and conjugating-transposing the system, we get 

Y^ = |X^D^|, (2) 

where Y = [yi,...,yp], X = [xi,...,xp] and denotes the conjugate-transpose of a ma¬ 
trix/vector. This reveals a “classic” phase retrieval problem: given the matrix of inputs X^, 
each column of Y^ is used to estimate each complex-valued column of D^. 

3.1. Phase retrieval 

The problem of reconstructing a complex vector given only the magnitude of measurements is 
a non-convex optimization problem notoriously difficult to solve. Many algorithms have been 
devised in the literature to deal with this problem. We can roughly divide them into three main 
families: 

1. The alternating-projection algorithms alternate projections on the span of the measure¬ 
ment matrix and on the object domain. Among these approaches, we can mention the 
works of Gerchberg & Saxton 1^ . Fienup ||28l and Griffin & Lim 1^ . 

2. The algorithms based on convex relaxations approximate the phase recovery problem by 
relaxed problems which can be solved efficiently by standard optimization procedures. 
Two of the main approaches of this type, namely PhaseLift 1^ and PhaseCut flM . rely 
in particular on semidefinite programming. 

3. The Bayesian approaches, recently envisaged in ||25]|26l, circumvent the non-linearity of 
the modulus through the introduction of hidden variables and resort to variational approx¬ 
imations to approximate the posterior distribution of the variables of interest. These latter 
methods have been shown to perform good reconstruction in a reasonable computational 
time 1261 . 

3.2. Bayesian variational approximations 

Additionally to the previous notations, we introduce new variables, modeling, on the one hand, 
the missing phases of the observations, and on the other hand, some acquisition noise. Thus, 
recalling that we resort to a conjugate-transposition of the matrix system, each absolute-valued 
measurement yjj, p G {1... F}, of any row y of Y, is expressed as 

+ (3) 

i=l 

where G [0,2;r) stands for its missing conjugate phase, x^i is the ith element of the pth row 
in X, J* corresponds to the ith conjugate element in the current estimated row d of D and COjx 
is an additive noise, assumed centered isotropic Gaussian (denoted the following) with 

variance G^. We moreover suppose that the probability distributions for the entries of the matrix 
and for the missing phases are: 


N 

i=\ 

and p{0) = Yi 
M=i 


(4) 


with 


(5) 


Under these assumptions, the absence of phases in the observations is naturally taken into 
account in the model since marginalizing on leads to a distribution on which only depends 
on the moduli of and T 4 L 1 d*. 

Within model ([^-([^, the recovery of the complex signal d can be expressed as the solution 
of the following marginalized Maximum A Posteriori (MAP) estimation problem 


d = argmaxp(d|y), 

d 

with p(d|y) = U(d,0|y). 

Jd 


( 6 ) 

(7) 


Because of the marginalization on the hidden variables 6 , the direct computation of /?(d|y) is 
however intractable in general. The solutions in |[25l[26]| optimally approximate, in a Kullback- 
Leibler sense, the posterior joint distribution p(d, 6 |y) by ^(d, 6 ) conditionally to a set of given 
constraints 

^(d,0) = argmin / A(d,0)log (8) 

JdJe Xd,0|y)' 

Depending on the minimization gives rise to different approximations. 

• In particular, ^ qi{di)Y[^=i defines a Mean-Field approxima¬ 

tion, and problem ^ can be efficiently solved using the “Variational Bayes Expectation- 
Maximization” (VBEM) algorithm ED . This is the approach considered in ll26l . denoted 
by prVBEM in the rest of this paper. 


With^ = {^|^ = 




} where [di... d^] (resp. [ 0 i ... 6 b ]) parti- 


tions the variables d (resp. 6 ), and at (resp. j3^) is the degree of variable node di (resp. 
0^), problem ^ refers to the minimization of the Bethe free energy, which can be solved 
by generalized approximate message passing (GAMP) algorithms, see (321. This is the 
approach followed in (251 . denoted by prGAMP in the rest of this paper. 


We will not detail here the structure of the resulting algorithms. We refer the interested reader 
to the papers (^ and (25i and the authors’ webpag^for a practical implementation of the 
prVBEM algorithm. 


3.3. Experiments and results 
3.3.1. Prediction performance 

To assess the accuracy of the TM estimated by the considered approaches, we adopt a cross¬ 
validation-like experimental framework. The setup is as follows. We measure the M = 40000 
camera pixels stemming from N = 900 DMD mirrors, 50% of them being turned on, the others 
off at each displayed pattern. The operation is repeated randomly P = 6000 times. Given this 
dataset, a row of the TM is then learned from p = aN calibration measurements, with a varying 
in {1,... ,6}, and used in a second step to predict the P — p remaining measurements. This 
estimation is performed on 50 different rows of the TM. 

We evaluate and compare the performance of 4 different algorithms: Gerchberg-Saxton (23, 
PhaseCut (24l, prGAMP (^ and prVBEM (26l. The algorithms present different complexi¬ 
ties. The implementation of PhaseCut (available on author’s webpag^ relies on interior-point 
methods, with a complexity growing as e)) where £ is the target precision (33l. 

Gerchberg-Saxton, prGAMP (in our own implementations) and prVBEM share similar com¬ 
plexities, of order ^{p^). 


^http://angelique.dremeau.free.fr/ (released October 7th, 2014) 
^http://www.di.ens.fr/data/software/ 
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Fig. 3: Prediction performance according to (a), the mean-square error (MSB), in log scale, 
and to (b), the normalized cross-correlation between observation predictions using the esti¬ 
mated TM, and actual measurements of the output moduli (square root of the camera intensity 
values), as a function of the number of calibration measurements (x-axis is a, such that p = aN 
calibration measurements are used). 


As a tradeoff between computational cost and performance, we set the stopping criteria 
for each algorithm as follows. PhaseCut is run until the target precision drops below 10“^. 
Gerchberg-Saxton stops after 3000 iterations, ;?rGAMP and prVBEM after 200 iterations. Given 
the complexities of the algorithms, we allow PhaseCut and Gerchberg-Saxon for a higher run¬ 
ning time, as shown and further discussed in Fig.[^ 

The prediction performance of the algorithms is evaluated according to the mean-square er¬ 
ror (MSB) (Fig. and the normalized cross-correlation (Fig. between the moduli of the 
P — p predicted measurements and the actual observed ones. For Gerchberg-Saxton, PhaseCut 
and prVBEM, the MSE curves present similar behaviors : they decrease monotonically with 
increasing a. This observation resonates in Fig.|^ with the general increase tendency of the 
correlation. Interestingly, we see that for a > 3, that is, for at least 3 times more real measure¬ 
ments than complex unknowns, prVBEM outperforms all other algorithms with a correlation 
around 0.95. 

On the contrary, prGAMP presents a contrasted performance. If it leads to a good corre¬ 
lation (Fig. between estimated and observed measurements, its MSE remains very high, 
independently of a (Fig. [^. Here, the algorithm finds an acceptable solution but only up to 
a multiplicative factor. While this is not an issue for the focusing experiment considered in 
Section]^ this could become a limitation in other more complex tasks. 

Finally, as previously mentioned, we allow in these experiments more iterations for PhaseCut 
and Gerchberg-Saxton. In parallel to the performance curves exposed above. Fig. [^illustrates 
the corresponding average running time of the 4 considered algorithms. In this figure, we see 
that prGAMP performs the lowest computational cost, closely followed by prVBEM. Phase- 
Cut requires a long running time, prohibitive in our application context. It should be noted 
that the Gerchberg-Saxton algorithm, although relatively slow, still exhibits good performance, 
especially given its simplicity. 

On the basis of these preliminary experiments, for the remaining of this paper, we choose 
the prVBEM approach and set the number of calibration measurements to p = AN, as a good 





















Fig. 4: Average running time (log-scale, in seconds) as a function of the number of calibra¬ 
tion measurements (x-axis is a, such that p = aN calibration measurements are used). These 
simulations have been done on a Macbook Air with a 1.7GHz i7 processor. 


tradeoff between performance and computation time. In this case, computing one row of the 
TM takes about .6 s in Matlab, on a Macbook Air with a 1.7GHz i7 processor, keeping in mind 
that rows are independent. 

3.3.2. Comparison of singular values to Random matrix theory 

Interestingly, we can check that the measured TM presents some characteristics as predicted 
by random matrix theory. One practical way is to verify that the distribution of its normalized 
singular values obeys the Marcenko-Pastur law 1341 . It should be noted that such apparently 
random signals are the hardest case for phase retrieval, where no specific structure can be taken 
into account. 

In order to reduce the infiuence of specifics of our experimental setting, we perform the 
following operations, as in |35l: 

i) We normalize over the rows and columns, to attenuate the illumination artifacts: residual 
illumination “by default” on each pixel of the camera for the rows, and inhomogeneous 
contribution of each DMD mirror on the entire set of camera pixels for the columns. 

ii) Because of the size of the speckle grains, two neighboring DMD mirrors may affect the 
material in the same way, as well, two pixels of the camera will be potentially correlated. 
To avoid this effect, we subsample the rows and columns of the matrix. 

To draw the empirical spectral density, we then consider the following setup. We subsample the 
columns of the matrix up to A = 200 and leave the number of rows varying, more precisely M = 
yN, with 7G {1,...,6}. These sub-matrices thus constitute partitions of the estimated matrix, 
randomly picked 100 times to average the resulting densities. Fig. [^compares the experimental 
curves to the theoretical ones drawn according to the Marcenko-Pastur law. We see that the 
experiments qualitatively follow the predictions. We remark, however that the larger y is, the 
more chances we have to consider the contributions of neighboring correlated pixels. This partly 
explains the increasing gap between both curves. 









Fig. 5: Density of the normalized singular values for different y = M/N. Stamped line: experi¬ 
mental results, continuous line: Marcenko-Pastur law. 


4. Focusing with the DMD 

Knowing the TM gives a powerful and flexible tool to control light within the scattering medium 
Ea . In particular, it can be used to compute which DMD input has to be set, in order to display 
a given arbitrary pattern at the receiver end. In this section, we demonstrate the special case of 
focusing light with maximum intensity on a desired pattern (a chosen sparse subset of the output 
pixels), with the TM measured experimentally as in the sectionIt should be emphasized that 
we keep the same experimental setup, with the binary DMD as input device. Here, simple 
inversion methods such as EH cannot be used, as these require intensity- or phase-modulated 
inputs. 

We propose here to resort to a similar Bayesian variational approach as for the calibration, 
adapted to the binary nature of the DMD inputs. 

4.1. Mean-Field-based inversion 

Formally, the problem can be expressed as an inverse problem, where, knowing the TM D and 
the observation y, we look for the DMD input x such as described in ([T]). Adopting a similar 
modeling as in previous section, we then assume, for all elements with /i G {1,... ,M}, 

N 

d^i Xi + Wjii ), (9) 

i=\ 

where G [0,2;r) stands for the missing conjugate phase, is the /ith element of the ith- 
column in D, X; G {0,1} corresponds to the state of the i-th DMD pixel and i s an additive 
noise, assumed centered isotropic Gaussian of variance cr^. As in the subsection |3.2[ we sup¬ 
pose that the elements 6 ^ are independently and uniformly distributed in the interval [0,2;r), 




























however, in order to accommodate for binary inputs, we consider here a Bernoulli model for x: 

P(x) = n^(^') p{xi)=Bcr{pi) = i^ txi = 0. 

Then, within model ([9l>-(fT0l>, we solve the marginalized MAP estimation: 

x = argmaxp(x|y), (11) 

X 

with p{x\y)= [ p{x,e\y), (12) 

Je 

and resort - following the comparison exposed in subsection |3.2| in the Gaussian case - to 
a Bayesian Mean-Field approximation. The particularization of the algorithm to the Bernoulli 
model is detailed in the appendix, an implementation is also available on author’s webpage. 

4.2. Experiments and results 



Fig. 6: Illustration of light focusing on 3 points. The circles mark the positions of the targets. 


In this section, we assess the performance of the proposed focusing approach through differ¬ 
ent experiments. The general setting is as follows. The DMD inputs, here taken of dimension 
N = 1600, are estimated according to the procedure described above from the desired outputs, 
focusing on 1 to 4 target points. We set the Bernoulli parameters p/ to 0.5, noticing that asymp¬ 
totically half of the DMD pixels are expected to be “ON” ( ifTTl ). Finally, the TM, reduced to its 
rows of interest, is measured as discussed in section 

Fig. shows an example of the observed output field, corresponding to the estimated DMD 
configuration, optimized to focus on 3 points. To quantitatively evaluate the focusing perfor¬ 
mance, we measure the intensity enhancement factor, as: 



(13) 


where /^c is the intensity inside the target area after spatial binary amplitude modulation is 
performed, /back is the average background intensity. This value is measured for 100 trials, as a 
function of the number of calibration measurements used to learn the TM. Two different setups 
are then considered: the single-point focusing case and the multi-target case. 






Fig. 7: Single target experiment. Enhancement factor as a function of the number of measure¬ 
ments used to learn the TM (x-axis is a, such that p = aN calibration measurements are used). 
For the same estimation of the TM, 2 focusing techniques are compared: phase conjugation 
(blue boxes), and the new Mean-Field technique (red boxes). 


4.3. Focusing on a single point 

Fig. [7] compares the enhancement factors achieved by two different focusing methods, namely 


a simple phase-conjugation - performing x = 


9t(D^y) > 0 


and the proposed method, in the 


case where only one target point is focused. Results are presented under a “box” format, where: 


• the middle segment stands for the average enhancement rj over the 100 trials, 

• the upper and lower bounds of the rectangle define the interval [f) — Orf ^FOrf] (where 

is the experimentally computed standard deviation), in which lies, under the Gaussian 
assumption, 68 % of the trials. 


• the whiskers represent the minimum and maximum values observed over the entire set of 
trials. 


For each experiment point a G {2,..., 6}, such that p = aN calibration measurements are used 
to compute the TM, we display the boxes related to the phase-conjugation method (blue boxes), 
and the new Mean-Field technique (red boxes) described in the previous section. 

As a first observation, we can see that the general dependency with regard to a noticeably 
resonates with the curve of the prVBEM algorithm in Fig.[^ there is a clear gap between the 
performance achieved for a = 2 and for a = 3, while, for a > 3, the intensity enhancement 
keeps increasing but less significantly. 

Interestingly, the Mean-Field approach seems to outperform phase-conjugation, with regard 
to the mean and maximum values measured, but not always in a statistically significant manner. 
Focusing on the most favorable case considered here, namely with a = 6, the best intensity 
enhancement factor lies around 140, to be compared with the ideal expected enhancement given 
by 1 + i (f - 1) - 255, see GTl- 

























4.4. Focusing on multiple points 

For this second setup, we are interested in the performance of the proposed algorithm in a 
context of multiple target points. Additionally to the intensity enhancement, we consider here 
the missed detection rate, defined as the number of trials (expressed in percentage) failing to 
focus on at least one of the multiple target points, i.e., the number of trials for which at least 
one of the T largest intensity peaks in the output image does not match any of the T targets. 



Fig. 8: Multiple target experiment, (a) Average enhancement factor as a function of the number 
of measurements used to learn the TM (x-axis is a, such that p = aN calibration measurements 
are used), and the number of target points (y-axis). (b) Missed detection rate (same axis as in 
(a)). 


Fig. [^represents these two figures of merit under diagram formats. They present an interest¬ 
ing general symmetry: increasing the number of targets or decreasing the number of calibration 
points leads to an increase of the missed detections and a decrease of the enhancement factor. 
The missed detection rate seems however more sensitive to the number of calibration points 
used to learn the TM: for a = 2 and 2 target points, the algorithm fails with a rate approach¬ 
ing 40%, while for a = 3 and the same number of targets, we keep a reasonable performance 
(around 10%). In a more general view, these figures greatly highlight the deep relation between 
the quality of the calibration and the focusing performance. 

5. Conclusion 

This paper shows that the full complex-valued transmission matrix of a strongly scattering mate¬ 
rial can be estimated, up to a global phase factor on each of its rows, with a simple experimental 
setup involving only real-valued inputs and outputs. In our experiment, the inputs are ampli¬ 
tude modulations on a binary DMD, and the output is the field intensity measured on a CCD 
camera, that gathers a significant amount of measurement noise. Note that no reference arm is 
used, that would allow interferometric measurements, but that would make the experimental 
setup more complex and considerably more unstable. 

We here resort to Bayesian phase retrieval techniques, and we have shown that, amongst such 
techniques, a recently proposed variational approach (VBEM) 1^ allows a precise estimation 
of the transmission matrix, tractable in computational complexity and scalable for large-size 
signals, provided that we have a sufficiently large number of input-output calibration signals. 











Experimental results validate this concept, both in terms of output prediction, distribution of 
singular values, and in an application of light focusing onto a number of target points in the 
output plane. It should be emphasized that this estimation of the transmission matrix opens 
many applications beyond light focusing, may it be for imaging through the scattering material 
1351 [36l, or for obtaining information about the scattering material itself. 

6. Appendix: Focusing with a Mean-Field based algorithm 

The VBEM algorithm is an iterative procedure which successively updates the factors of the 
Mean-Eield approximation. Particularized to model (|9|)-(p^, this gives raise to the following 
update equations: 






exp 




^{xi) = p{xi) exp ( X, 


2 9t(df{r,))-dfd, 


where 


(ri) =y-E^(^'t= l)d;t 

kjti 




y = 

= Y^^ixi = 1 ) d^i, 


(14) 

(15) 

(16) 

(17) 

(18) 


and lo (resp. Ii) stands for the modified Bessel function of the first kind for order 0 (resp. 1). 
Coming back to problem ( pTj ), an approximation of p(x|y) thus simply follows from 


p{My) = f p(x,0|y), 

Jd 

= Y[q{xi). 


(19) 

( 20 ) 
( 21 ) 


Using this approximation, the problem is then easy to solve by a simple thresholding operation, 
i.e., Xi = lif q{xi = 1) >0.5 and Xi = 0 otherwise. 
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